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Abstract 

The quadratic contact process (QCP) is a natural extension of the well studied linear contact 
process where infected (1) individuals infect susceptible (0) neighbors at rate A and infected in- 
dividuals recover (1 — > 0) at rate 1. In the QCP, a combination of two I's is required to effect 
a — > 1 change. We extend the study of the QCP, which so far has been limited to lattices, to 
complex networks. We define two versions of the QCP - vertex centered (VQCP) and edge centered 
(EQCP) with birth events 1 — — 1 — > 1 — 1 — 1 and 1-1-0 — > 1 — 1 — 1 respectively, where 
'— ' represents an edge. We investigate the effects of network topology by considering the QCP on 
random regular, Erdos-Renyi and power law random graphs. We perform mean field calculations 
as well as simulations to find the steady state fraction of occupied vertices as a function of the birth 
rate. We find that on the random regular and Erdos-Renyi graphs, there is a discontinuous phase 
transition with a region of bistability, whereas on the heavy tailed power law graph, the transition 
is continuous. The critical birth rate is found to be positive in the former but zero in the latter. 
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I. INTRODUCTION 



Inspired by technological and social networks, the study of complex networks has seen a 
surge in the past fifteen years |lH5]. Research has traditionally progressed in two distinct 
directions - dynamics of networks and dynamics on networks. The former is concerned with 
the formation of a network or change in its structure with time, whereas the latter deals 
with processes (deterministic or stochastic) taking place on a fixed network. Preferential 
attachment and its many generalizations |6l |7] are prototypical examples of the first type. 
Examples of the second are epidemics [HHTT], the voter model for the spread of an opinion 
[T2HT^ . cascades [T5HT7] that model spread of a technology, and evolutionary games [T8] . 
The phase transitions [IHl [20] associated with these models have been of particular interest. 

In the mathematics community, spatial models are studied under the heading of interact- 
ing particle systems [21] . One of the simplest models of those models is the contact process 
[22H21] (equivalent to the SIS model in epidemiology). In the hnear contact process each site 
can be in one of two states which we will call 1 and 0. O's become 1 at a rate proportional 
to the number of 1 neighbors they have and I's become at a constant rate (here and in all 
following models, unless otherwise specified, the processes occur in continuous time). 

A natural extension of the linear process is the quadratic contact process (QCP) where 
each — 1 event will require two other sites in state 1. We will occasionally refer to 1 as 
being the "occupied" state and as being "vacant", and the events — )■ 1 and 1 — t- to 
be birth and death events respectively. At this stage, the model is quite general in that we 
do not specify where the two I's that cause the — )■ 1 event must be located with respect 
to the 0. On the 2D lattice, specifying these locations leads to different realizations of the 
QCP. For example, Toom's North-East-Center model (originally defined in discrete time) 
allows a at site x to be filled if its neighbors x + (0, 1) and x + (1, 0) are occupied |25] . 
Chen [261 122] has studied versions of Toom's model in which two or three specified adjacent 
pairs or all four adjacent pairs are allowed to reproduce. Evans, Guo and Liu [2HHS2] have 
studied the QCP as a model for adsorption-desorption on a two dimensional square lattice. 
In the version of the model studied by Liu [32], becomes 1 at rate proportion to the 
number of adjacent pairs of 1 neighbors. He found a discontinuous phase transition with 
a region of bistability, where the I's die out starting from a small density. He also found 
that by introducing spontaneous births at a sufficiently high rate, the transition becomes 
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continuous. 

The QCP is similar to Schlogl's second model of autocatalysis characterized by 
chemical reactions 2X — > 3X,X — > where X represents the reactant. Grassberger [31] 
studied a version of Schlogl's second model in which each site has a maximum occupancy 
of two and doubly occupied sites give birth to a neighboring vacant site. He found that the 
model shows a continuous phase transition in 2D. 

Studies to date on the QCP have been limited to regular lattices in low dimensions. In 
this paper, we extend the study to complex networks. There are two ways to view the QCP 
on networks: 

• as a model that replaces the linear birth rate of the contact process that has been 
extensively studied on networks [35l |36], by a quadratic birth rate. 



• as an alternative model for the spread of rumors, fads and technologies such as smart 
phones in a social network. In sociology the requirement of more than a single 1 for the 
"birth" event is called complex contagion |37]. Also related are the threshold contact 
process |3H] and models for the study of "cascades" jTB]. The key difference here is 
that the QCP involves a death event that represents the loss of interest in the fad or 
technology and the rate for birth events is a function of the actual number and not 
the fraction of occupied neighbors. 

The questions we are interested in are: How does network topology affect these phase 
transitions? What model and network features lead to discontinuous versus continuous 
phase transitions? 

The paper is organized as follows. We define the specific QCP that we study in Section 



O and we do mean field calculations in Section HI In Section IV we present a few rigorous 



results about the QCP. Simulation results are presented in Section IVl followed by some 



concluding remarks in Section VI 



II. MODEL DEFINITION 



The birth event in the linear contact process can be formulated as each 1 — edge converts 
to a 1 — 1 edge at a constant rate A. Such a definition can be easily extended to the quadratic 
case by defining the birth event in terms of connected vertex triples. Two such definitions 
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are possible: 1 — — 1 — > 1 — 1 — 1, and 1-1-0 — > 1 — 1 — 1. We call the former version 
the vertex centered QCP (VQCP) because the central vertex is getting filled by its two 
neighboring Is, and the latter as the edge centered QCP (EQCP) as it can be viewed as a 
1 — 1 edge giving birth on to a neighboring vacant vertex. Note that the models can also 
be defined in terms of how a vacant vertex gets filled i.e., suppose that a vertex has k 1 
neighbors and j 1 — 1 neighbors [39], then the vertex will become 1 at rates (2) A and jX 
in the VQCP and EQCP respectively. Death events 1 — > occur dbt Tcttc 1 as in the linear 
process. 

If the death rate is changed to zero, the VQCP reduces to bootstrap percolation [ID] 
where vertices that are occupied remain occupied forever and vacant vertices that have at 
least two occupied neighbors become occupied. While bootstrap percolation is typically 
defined in discrete time, the final configuration of the network is independent of whether the 
dynamics happens in discrete or continuous (as in our model) time. 

We will use random graphs as models for complex networks on which the QCP is taking 
place. We will denote by d the degree of a randomly chosen vertex in the network and the 
degree distribution by pk = P(c? = k). We are interested in networks with size n — )■ 00 and 
where the vertex degrees are uncorrelated. The specific random graphs that we will consider 
are 

• Random regular graphs RR(/i) in which each vertex has degree fi. Since everyone has 
exactly /i friends, this graph is not a good model of a social network. However, the 
fact that it looks locally like a tree will facilitate proving results. 

• Erdos-Renyi random graphs ER(/i) where each pair of vertices is connected with prob- 
ability /i/n. In the n — )■ 00 limit, the degree distribution of the limiting graph is 
Poisson with mean /x. This is a prototypical model for the situation in which the 
degree distribution has a rapidly decaying tail. 

• Power law random graphs PL (a) with degree distribution p^ = ck~°'. We are particu- 
larly interested in graphs where the exponent a lies between 2 and 3, which has been 
found to be the case for many real world networks |41j . 

We will occasionally refer to RR and ER as homogeneous networks as their degree distri- 
butions are peaked around the mean, in contrast to PL where the distribution has a heavy 
tail. 
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III. MEAN FIELD CALCULATIONS 



We can attempt an analytical study of the dynamics by writing the equations for the 
various moments of the network. Let g he a small graph labeled with I's and O's. We define 
the (7— moment, written as (g), of a {0,1} valued process on a graph G as the expected 
number of copies of g that exist in the set of all subgraphs of G. For example if g = 1 — — 1, 
we look at all the connected vertex triples in the network and count the ones where the center 
vertex is in state and the other two vertices are in state 1. We will write p(X,p^^^;t) as 
the density (1)/"^. at time t with a birth rate of A and an initial configuration where each 
vertex is independently occupied with a probability p^^\ The order parameter for our phase 
transitions is the steady state density 

p.(A,p(°))= limp(A,p(°);t). (1) 

We define the critical birth rate Ac as the birthrate above which there exists a stable steady 
state density that is greater than zero, i.e., 

Ac = inf{A:p.(A,l) >0}. (2) 

In the definition above, we chose p^"-* = 1 since it has the best chance of having a positive 
limit. We also define a critical initial density pc as the minimum initial density required to 
reach a positive steady state density when the birth rate is infinite, i.e., 

Pc = inf{p(°): limp,(A,p(°))>0}. (3) 

A— >CXD 

From their definitions, it is straight forward to write the dynamical equations of (1) for 
the VQCP and the EQCP, 

d I (1 - - 1) for the VQCP 

;n(l) = -(l) + Ar . (4) 

1(1-1-0) for the EQCP 

If we were to write the equations for the third order moments that appear on the RHS of 
Q, those equation would involve still higher order moments. Continuing this way, we end 
up with an infinite series of equations that are not closed. Therefore we resort to a mean 
field approximation by assuming the states of neighbors of a vertex to be independent at all 
times. 
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FIG. 1. The blue, dashed and red curves correspond to p,, = 0, />_ and p+ respectively obtained 
from the mean field calculation for both QCP types on homogeneous networks. 

A. Homogeneous networks 



In the following we do a naive calculation that ignores the correlation between degree and 
occupancy, which should be reasonable for homogeneous networks. With these assumptions, 
(1 — — 1) will be np^(l — p)((2))- Plugging this value into Q we get 



P 



(5) 



Setting the RHS of ([s]) to zero gives a cubic equation whose roots are the possible steady 
state densities p^. Clearly, zero is a trivial root of (|5]). The other two roots are 



P± 



1±A 1 



A 



(6) 



These solutions are real only when A > Ac = 4/ ((2))- In the language of nonlinear dynamics, 
(|5| exhibits a saddle node bifurcation at Ac. It is easy to see that zero and p^ are stable 
fixed points whereas p_ is an unstable fixed point. This can be seen in Fig. [TJ The limiting 
critical initial density is 

Pc = lim p_ = . (7) 

A— >oo 

For ER(p) we have ((2)) = which gives Ac = For PL(a; < 3) we have 

((2)) = cxD so Ac = while PL(a > 3) has finite ((2)) leading to a non-zero value for Ac. 
The mean field calculation for the EQCP is essentially the same as done above and predicts 



the same qualitative features. Thus, for networks with finite ((2)), the simple mean field 
calculation predicts a discontinuous phase transition at A = Ac and a region of bistability 
for A > Ac, for both QCP types. 



B. Heavy tailed degree distributions 



The mean field calculation of Section III A is simplistic since it ignores the fact that the 
occupancy probability depends on the degree. Pastor-Satorras and Vespignani |H] improved 
the mean field approach for the linear contact process by defining p^, the fraction of vertices 
of degree k that are occupied, and 9, the probability that a given edge points to an occupied 
vertex. These variables can be related through the size biased degree distribution qk = 
kpk/ (d) which is the distribution of the degree of a vertex at the end of a randomly chosen 
edge. 

= J2lkPk (8) 

k 

Note that for homogeneous networks we assumed 6 = p. As before, the state of the neighbors 
of a vacant vertex are assumed to be independent. So the number of occupied neighbors 
of a vertex of degree k follow a distribution Binomial(/c, ^). This enables us to apply this 
approach to the VQCP. We write equations for p^, 

'k 



Pk 



-Pfc + A(l-pfe)( 2 \9' 



So in steady state 



Pk* 



1 + A©02 



(9) 



(10) 



Combining m\ and (10) leads us to a self-consistent equation for 9^,. 



where 



/(A, 



9* = 9^:I{\, 9^) 



kpk 
{d) 



XQ9 



\Q9^ 



(11) 



(12) 



Clearly, = is a solution of (11). Finding a non-trivial solution involves solving 



J(A,^,) = 1 , G(0,1). 



(13) 



For power law graphs PL (a), as we will now show, the mean field calculation predicts 
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FIG. 2. /(A, 6) versus 6 near A = Ac for various power law graphs. 



If a > 3, Ac > and the transition is discontinuous. 



If a = 3, Ac > and the transition is continuous. 



• If 2 < a < 3, Ac = and the transition is continuous. 
To begin we note that if pk ~ Ck^°' with a > 4 then 

ex 

W) 



k 



as ^ — 0. When 3 < a < 4 breaking the sum at A; = 1/^, and bounding the fraction hj C/9 
for k >l/6 (here and in what follows C is some finite constant). 

as — 7- 0. Since J(A, 1) < 1 for all A, the curve will have supg^jg,!] -^(Ac, ^) = 1 at some Ac. 
Then for A > Ac there will be two roots, the larger of which is the relevant solution. 
When 2 < a < 3 approximating the sum by an integral we have 

/(A,<))~r-^^-^d.. (16) 

From this we see that if 2 < a < 3 then as ^ — )• 0, I{X,6) — )■ oo, so there is a solution to 
I{X,dx) for any A > 0, and — as A — ?■ 0. For small values of p^:, a little calculation 
shows that the following relation holds 

~ (A - Ac)^(") where j{a) = - ^ (17) 

6 — a 2 



is the critical exponent. 

In the borderhne case a = 3 the hmit as — t- is finite. There is a critical value Ac so 
that /(Ac, 0) = 1. Since 9 — > /(A, 9) is decreasing, it follows that for A > Ac we have one 
solution /(A, 9\) = 1 which has — as A — ?■ Ac- The behavior of /(A ~ Ac, 9) is shown in 
Fig.|2} 

A second way to determine the nature of the phase transition is to adapt the argument 
of Gleeson and Cahalane (TU], which can be applied if we use a discrete time version of the 
model in which a vertex with k neighboring pairs will be occupied at the next step with 
probability l — {l—p)^. The computation in their formulas (l)-(3) supposes that the vertices 
at a distance n from x are independently occupied with probability po- The function G{p) 
defined in their (3) gives the occupancy probabilities at distance k — 1 assuming that the 
probabilities at distance k are p. Iterating G n times and letting n ^ oo gives a prediction 
about the limiting density in the cascade, if one repeats the calculation for our system then 
is an unstable fixed point when a < 3, while it is locally attracting for a > 3. This agrees 
with the mean-field prediction of Ac = in the former case and a discontinuous transition 
with Ac > in the second. 

IV. SOME RIGOROUS RESULTS 

We have not been able to extend the mean field calculation to the EQCP on power law 
graphs, but by generalizing an argument of Chatterjee and Durrett [35j we can prove that 
Ac = for a G (2,oo). The details are somewhat lengthy, so we only explain the main 
idea. Consider a tree in which the vertex has k neighbors and each of its neighbors has 
/ neighbors and / is chosen so that l\ > 10. One can show that if k is large then with 
high probability the infection will persist on this graph for time > exp(c(A)A;). In a power 
law graph one can find such trees with k = n^/^°'~^\ Using the prolonged persistence on 
these trees as a building block one can easily show that if we start with all vertices occupied 
the infection persists for time > exp(r;,^~'^) with a positive fraction of the vertices occupied. 
With more work (see ^2] B3] ) on can prove persistence for time exp(c(A)n). 

For both types of QCP it is easy to show that it is impossible to have a discontinuous 
transition with Ac = 0. The proof for VQCP is as follows. Let (Ik) be the expected number 
of occupied sites of degree k and (lOfcl) be the expected number of 1-0-1 triples when the 
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vertex has degree k. We can write an equation similar to Q 

d 



^^.Ifc) = -(lfc) + A(l-Ofc-l) (18) 



which means at steady state 

(Ifc), = A(l-Ofc-l), <A(Ofc). 



Pk* < A(l - pk*) { ^] ^ Pk* < -T^r • (19) 



1 + A© 



So, as A — )■ 0, Pk* and p* = J2k Pk*Pk — ^ 0. Thus the transition will be continuous. The 
proof for EQCP is similar. In that case the subscript k stands for the secondary degree 
which is defined as the number of neighbors of neighbors of a given vertex (not including 
itself), i.e., S'^\x) = \{z : z y,y x, z ^ x}\. 

^(U) = -(lfc) + A(l-l-Ofc) (20) 



So at steady state 



Xk 

iU)* = A(l - 1 - Ofc), < X{Ok)*k < (21) 



Thus for both QCP types we find that if Ac = then the phase transition is continuous. 

For both QCP types on random r-regular graphs we can show that the critical birth 
rate is positive as follows. In the EQCP let there be m occupied vertices. Each of these m 
vertices can have at most r neighbors that are vacant and can give birth on to them at rate 
< (r — 1)A or die at rate 1. So the total birth rate in the network is < (r — l)Arm against a 
death rate of m, and it follows that Ac > l/r(r — 1). Similarly, for the VQCP the total birth 
rate is < X(^Qrm and it follows that Ac > l/r(^Q. These arguments depend on the degree 
being bounded, so they do not work for Erdos-Renyi and power law graphs. 



V. SIMULATION RESULTS 



We perform simulations of the QCP on RR(4), ER(4) and PL (2. 5). We generate the 
random regular and power law random graphs using the recipe called configuration model 
|44j . We draw samples dx from the degree distribution and attach that many "half-edges" to 
vertex x. We pair all the half edges in the network at random. We then delete all self loops 
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FIG. 3. Steady state density reached, starting from all vertices occupied, for QCP on homogeneous 
networks of various sizes n. 

and multiple edges. When a > 2 this does not significantly modify the degree distribution. 
If turns out to be odd (an event with probability ^ we ignore the last remaining 

unpaired half-edge. Furthermore, for PL (2. 5) we start the degree distribution at 3 as, in the 
VQCP, the vertices of degree 1 and 2 are impossible or difficult to get occupied. 

To deal with finite size effects, we observe how the plot of the steady state density p*(A, 1) 
versus A starting with all vertices occupied changes when size n of the network ranging from 
10^ to 10^. Fig. [sjshows the results of both QCP types on RR(4) and ER(4). Here the curves 
seem to converge to a positive value implying a positive Ac. The results for PL (2. 5) are shown 
in Fig. |4| We observe that the transition happens close to zero and moves towards zero with 
increasing n indicating that the critical birth rate is zero. As explained earlier, if Ac = 
then the transition is continuous. This is consistent with the the mean field predictions for 
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the VQCP and rigorous result for EQCP. In addition, in Fig. |4[a), the critical exponent for 
the n = 10^ curve can be measured to be approximately 1.45 which is close to the mean 
field value of 1.5 (obtained by setting a = 2.5 in ([It))). 

In order to further investigate the phase transitions in random regular and Erdos-Renyi 
graphs, we look at the steady state density attained by starting from two different initial 
densities for the same network size n = 10^. Fig. [s] again shows a similar pattern across 
both QCP and both network types. In Fig. ^h) we see that for birth rates between 0.9 
and 2.3, the VQCP survives when the starting configuration had all vertices occupied but 
dies out when starting with only one-tenth of the vertices occupied. Thus we see bistability 
in the region A G (0.9, 2.3) implying a discontinuous transition and consequently that Ac is 
positive and close to 0.9. This is qualitatively in agreement with the mean field prediction 
seen in Fig. [T| although the critical birth rate of 0.9 shows a deviation from the mean field 
value of 8/42 = 0.5. 

Pontes and Schonmann [IS] have shown that for bootstrap percolation on the tree there is 
a critical density pc so that if the initial density is < pc then the final bootstrap percolation 
configuration has no giant component of occupied sites. In this situation having deaths at a 
positive rate in the VQCP will lead to an empty configuration. The last argument is for the 
tree, but results of Balogh and Pittel jl6] show that similar conclusions hold on the random 
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(d) EQCP on ER(4) 
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FIG. 5. Steady state density reached, starting from two different initial densities p^'^^ for QCP on 
homogeneous networks of size n = 10^. Notice the similarity with the mean field prediction on 
Fig.[Tj 



regular graph. While this argument is not completely rigorous, the reader should note that 
since all of the VQCP are dominated by bootstrap percolation, it follows that the limiting 
critical initial density defined in ([s]) has pc > in contrast to the mean field prediction in 
(j?]). Fig. [6] shows the final density attained as a function of the initial density when the 
birth rate is infinite (and death rate is positive). We see that pc in the VQCP is positive 
for the random regular and Erdos-Renyi graphs whereas it is zero for the power law graph. 
The corresponding results (not shown here) in the case of the EQCP indicate that pc = 
for random regular, Erdos-Renyi and power law random graphs. 
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FIG. 6. Steady state density when the birth rate is infinite in the VQCP. Note that the p*-"-* axis 
is in the log scale. 
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cont. , + HZlllSj 


NA , oo 


cont. , + [19] 


2D 


cont. , + 


discon. , + 


cont. , + [S! 


RR 


cont. la] , + [52] 


discon. , + 


discon. , + 


ER 


cont. , + [36j 


discon. , + 


discon. , + 


PL((2,3)) 


cont. , |35] 


cont. , 


cont. , 


PL(3) 


cont. , ^ 


cont. , + 


cont. , 


PL((3,oo)) 


cont. , [35] 


discon. , + 


cont. , 



TABLE L Nature of phase transitions of contact processes on various networks. Note that '0','+' 
and 'cx)' stand for zero, positive and infinite values respectively of Ac. 

VI. CONCLUSION 

In this paper we have investigated the properties of two versions of the quadratic contact 
process on three types of random graphs. The mean field calculations we performed agree 
with the simulation results. This may be due to the fact that complex networks have 
exponential volume growth and therefore are like infinite dimensional lattices where mean 
field is exact. 

Table |I] summarizes what is known about the phase transitions of contact processes in 
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1 and 2 dimensional lattices and on the random graphs RR, ER and PL. The positivity of 
the critical birth rate for ID, 2D and RR follows trivially from the boundedness of their 
degrees. For VQCP on a ID lattice, two consecutive O's can never get filled and it follows 
that Ac = oo. The results for the linear process on RR are inferred from the rigorous results 
for trees and the fact that RR is locally tree like. 

The results indicate that the EQCP is qualitatively not very different from the linear 
contact process on low dimensional lattices and power law graphs, in contrast to the VQCP, 
which differs from its low dimensional analogue. In view of the fact that they are very 
different in how they fill vacant vertices on a network, the similarity between VQCP and 
EQCP in their phase transitions on complex networks is a little perplexing. 

The EQCP can easily propagate on a chain and "cross bridges" connecting communities, 
compared to the VQCP which always requires two occupied neighbors. In the EQCP vertices 
with a large number of neighbors of large degree are the key to its survival. However, in the 
VQCP it is impossible for the central vertices to repopulate the leaves, so these structures 
are not long lasting. In contrast the Gleeson-Cahalane calculation suggests that survival is 
due to the fact that as waves of particles move through the system the densities increase. 

It is also interesting to note that the discontinuous phase transitions in the QCP are ac- 
companied by regions of bistability. Such a feature seems to be typical of all non-equilibrium 
phase transitions [531 El] • 
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